function [shortTermSeq, repeat_yesno, guideFunction] = ProbabilityRepetition_sequence(seqLen,period_labels,plot_boolean)

if nargin < 3
    plot_boolean = 0;
end

plot_boolean = logical(plot_boolean);

if nargin < 2
    period_labels = 7:15;
end

if nargin < 1
    seqLen = 128;
end


if mod(seqLen,2) ~= 0
    error('seqLen should be even');
end


% Generate guide function
guideFunction = zeros(1,seqLen);

for thisPeriod = period_labels
    guideFunction = guideFunction + sin((2*pi/thisPeriod)*(1:seqLen) + (rand * pi));
end

guideFunction = guideFunction - min(guideFunction);
guideFunction = guideFunction ./ max(guideFunction);


% Generate sequence of repeats
repeat_yesno = zeros(1,seqLen);

for i=1:seqLen
    repeat_yesno(i) = binornd(1,guideFunction(i));
end

% Generate the real sequence
shortTermSeq = zeros(1,seqLen);
% Randomize the first element (doesn't matter, since it' symmetric)
shortTermSeq(1) = binornd(1,0.5);
for i=2:seqLen
    if repeat_yesno(i) == 1
        shortTermSeq(i) = shortTermSeq(i-1);
    else
        shortTermSeq(i) = mod(shortTermSeq(i-1)+1,2);
    end
end

% Adjust the first element of repeat_yesno
if shortTermSeq(1) == shortTermSeq(end)
    repeat_yesno(1) = 1;
else
    repeat_yesno(1) = 0;
end


% Adjust the sequence so that it has the same number of 1s and 0s
excessOfOnes = sum(shortTermSeq) - (seqLen/2);
[~,sequenceOfIndicesToChange] = sort(abs(guideFunction-0.5),'descend');

index_to_change = 0;
while excessOfOnes ~= 0 && index_to_change < seqLen
    
    index_to_change = 1;
    % If there are more 1s than 0s, remove some 1s
    if excessOfOnes > 0
        % Find the preferred index to flip
        while (repeat_yesno(sequenceOfIndicesToChange(index_to_change)) == round(guideFunction(sequenceOfIndicesToChange(index_to_change)))) || (shortTermSeq(sequenceOfIndicesToChange(index_to_change)) == 0)
            index_to_change = index_to_change + 1;
            if index_to_change == seqLen
                break
            end
        end
        % Flip it!
        shortTermSeq(sequenceOfIndicesToChange(index_to_change)) = 0;
        
    else % If there are more 0s than 1s, remove some 0s
        % Find the preferred index to flip
        while (repeat_yesno(sequenceOfIndicesToChange(index_to_change)) == round(guideFunction(sequenceOfIndicesToChange(index_to_change)))) || (shortTermSeq(sequenceOfIndicesToChange(index_to_change)) == 1)
            index_to_change = index_to_change + 1;
            if index_to_change == seqLen
                break
            end
        end
        % Flip it!
        shortTermSeq(sequenceOfIndicesToChange(index_to_change)) = 1;
        
    end
    
    % Re-generate repeat_yesno sequence
    for i=2:seqLen
        if shortTermSeq(i) == shortTermSeq(i-1)
            repeat_yesno(i) = 1;
        else
            repeat_yesno(i) = 0;
        end
    end
    if shortTermSeq(1) == shortTermSeq(end)
        repeat_yesno(1) = 1;
    else
        repeat_yesno(1) = 0;
    end
    
    % Calculate the new excess of ones
    excessOfOnes = sum(shortTermSeq) - (seqLen/2);
end


%{
if sum(shortTermSeq==1) ~= sum(shortTermSeq==0)
    error('Something went wrong. Check the code');
end
%}

if plot_boolean
    plot(guideFunction)
    hold on
    plot(shortTermSeq,'r')
end





